!> author: 左志华
!> date: 2022-07-05
!>
!> Artificial viscosity and heat <br>
!> 人工热量与人工粘性
!> @todo 人工热量，人工热量影响内能
module sph_artifical_force

    use sph_kinds, only: rk
    use sph_particle, only: particle_t
    use twin_sph_list_type, only: twin_list_t
    use sph_env, only: nums

    logical :: art_visc                     !! Artificial viscosity <br>
                                            !! 人工粘度
    logical :: art_heat                     !! Artificial heat <br>
                                            !! 人工热量

    real(rk), parameter :: alpha = 1.0_rk   !! 剪切粘度
    real(rk), parameter :: beta = 1.0_rk    !! 容积（bulk）粘度
    real(rk), parameter :: etq = 0.1_rk     !! 避免奇点

contains

    !> Artificial viscosity <br>
    !> 人工粘度
    !> @todo 添加测试
    pure subroutine artificial_viscosity(particle, pairs, acc)
        type(particle_t), intent(inout) :: particle     !! particle <br>
                                                    !! 计算域
        type(twin_list_t), intent(in) :: pairs      !! Pair list <br>
                                                    !! 粒子对
        real(rk), intent(inout) :: acc(:, :)        !! Visc acceleration <br>
                                                    !! 粘性力（加速度）
        real(rk) :: vr, rr, dvel(2), dx(2), h(2)
        real(rk) :: muv, mc, mrho, piv
        integer :: i

        do concurrent(i=1:pairs%log)
            associate (ip => pairs%items(1, i), jp => pairs%items(2, i))

                if (ip > nums(1) .or. jp > nums(1)) cycle         ! 虚粒子不参与内力计算

                dvel(:) = particle%vel(:, ip) - particle%vel(:, jp)
                dx(:) = particle%loc(:, ip) - particle%loc(:, jp)
                vr = sum(dvel(:)*dx(:))
                rr = sum(dx(:)**2)

                if (vr < 0) then
                    muv = particle%hsml*vr/(rr + (particle%hsml*etq)**2)
                    mc = (particle%c(ip) + particle%c(jp))/2
                    mrho = (particle%rho(ip) + particle%rho(jp))/2
                    piv = (beta*muv - alpha*mc)*muv/mrho

                    h(:) = -piv*pairs%dwdx(:, i)
                    acc(:, ip) = acc(:, ip) + h(:)*particle%mass(jp)
                    acc(:, jp) = acc(:, jp) - h(:)*particle%mass(ip)
                end if

            end associate
        end do

    end subroutine artificial_viscosity

end module sph_artifical_force
